                                                                       /* In His Name */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>

#ifndef SIZE
#define SIZE 1024 //size of rows.
#endif

#ifndef THRESHOLD
#define THRESHOLD 64 
#endif

/** 
 * Assum all matrices are square and the number of rows is 
 * aslways a multiple of 2.
 **/

// Add two matrices
void AddMatrix(float* a, float* b, float* c, int size);
// Subtrace Two Matrices
void SubMatrix(float* a, float* b, float* c, int size);
// Naively Multiple Two Matrices 
void NaiveMul(float* a, float* b, float* c, int size);
// Strassen's algorithm
void Strassen(float* a, float* b, float* c, int size);
// Print a matrix
void printMatrix(float* a, int size);

int main(int argc, char** argv) {
  float* a, * b, *c;
  a = (float*)(malloc(SIZE*SIZE*sizeof(float)));
  b = (float*)(malloc(SIZE*SIZE*sizeof(float)));
  c = (float*)(malloc(SIZE*SIZE*sizeof(float)));

  int i = 0 , j = 0;
  for(; i < SIZE ; i++) {
    j = 0;
    for(; j < SIZE ; j++) {
      a[i*64 + j] = 1;
      b[i*64 + j] = 3;
    }
  }
  memset(c,0,SIZE*SIZE*sizeof(float));

  printf("---> Using the Naive implementation!\n");
  timeval tim;
  gettimeofday(&tim, NULL);
  double t1=tim.tv_sec+(tim.tv_usec/1000000.0);	
  NaiveMul(a,b,c,SIZE);
  //printf("A is \n");
  //printMatrix(a,64);
  //printf("----------------------------------\n");
  //printf("B is \n");
  //printMatrix(b,64);
  //printf("----------------------------------\n");
  //printf("A*B is \n");
  //printMatrix(c,64);
  gettimeofday(&tim, NULL);
  double t2=tim.tv_sec+(tim.tv_usec/1000000.0);	
  printf("%.6lf seconds elapsed\n", t2-t1);
  memset(c,0,SIZE*SIZE*sizeof(float));

  printf("---> Using Strassen's Algorithm\n");
  gettimeofday(&tim, NULL);
  t1=tim.tv_sec+(tim.tv_usec/1000000.0);	
  Strassen(a,b,c,SIZE);
  gettimeofday(&tim, NULL);
  t2=tim.tv_sec+(tim.tv_usec/1000000.0);
  printf("%.6lf seconds elapsed\n", t2-t1);	

  free(a);
  free(b);
  free(c);
  
  return 1;
}

void AddMatrix(float* a, float* b, float* c, int size) {
  int i,j;
  for(i = 0 ; i < size ; i++) {
    for (j = 0 ; j < size ; j++) {
      c[i*size + j] = a[i*size + j] + b[i*size + j];
    }
  }
}

void SubMatrix(float* a, float* b, float* c, int size) {
  int i,j;
  for(i = 0 ; i < size ; i++) {
    for (j = 0 ; j < size ; j++) {
      c[i*size + j] = a[i*size + j] - b[i*size + j];
    }
  }
}

void NaiveMul(float* a, float* b, float* c, int size) {
  int i, j, k;
  size_t _size = size*size*sizeof(float);
  memset(c,0,_size);
  for(i = 0 ; i < size ; i++) {
    for(k = 0 ; k < size ; k++) {
      for(j = 0 ; j < size ; j++) {
	c[i*size + j] += a[i*size+j]*b[k*size + j];
      }
    }
  }
}

void Strassen(float* a, float* b, float* c, int size) {
  if (size <= THRESHOLD) {
    return NaiveMul(a,b,c,size);
  }

  int half = size/2;
  size_t _size = half*half*sizeof(float);

  /* Allocate the need submatrices */
  float* a11,* a12,* a21,* a22,
    * b11,* b12,* b21,* b22,
    * c11,* c12,* c21,* c22,
    * M1,* M2,* M3,* M4,* M5,* M6,* M7,
    * atemp, * btemp;

  a11 = (float*)(malloc(_size));
  a12 = (float*)(malloc(_size));
  a21 = (float*)(malloc(_size));
  a22 = (float*)(malloc(_size));
  b11 = (float*)(malloc(_size));
  b12 = (float*)(malloc(_size));
  b21 = (float*)(malloc(_size));
  b22 = (float*)(malloc(_size));
  c11 = (float*)(malloc(_size));
  c12 = (float*)(malloc(_size));
  c21 = (float*)(malloc(_size));
  c22 = (float*)(malloc(_size));
  M1 = (float*)(malloc(_size));
  M2 = (float*)(malloc(_size));
  M3 = (float*)(malloc(_size));
  M4 = (float*)(malloc(_size));
  M5 = (float*)(malloc(_size));
  M6 = (float*)(malloc(_size));
  M7 = (float*)(malloc(_size));
  atemp = (float*)(malloc(_size));
  btemp = (float*)(malloc(_size));

  /* Fill the matrices appropriately */
  int i , j;
  for(i=0 ; i < half ; i++) {
    for(j = 0 ; j < half ; j++) {
      a11[i*half + j] = a[i*size + j];
      a12[i*half + j] = a[(half+i)*size + j];
      a21[i*half + j] = a[i*size + (half+j)];
      a22[i*half + j] = a[(half+i)*size + (half+j)];
      b11[i*half + j] = b[i*size + j];
      b12[i*half + j] = b[(half+i)*size + j];
      b21[i*half + j] = b[i*size + (half+j)];
      b22[i*half + j] = b[(half+i)*size + (half+j)];
    }
  }

  /* Start Applying the algorithm */

  // M1 = (A11 + A22)(B11 + B22)
  AddMatrix(a11,a22,atemp,half);
  AddMatrix(b11,b22,btemp,half);
  Strassen(atemp, btemp, M1, half);

  // M2 = (A21 + A22)B11
  AddMatrix(a21, a22, atemp, half);
  Strassen(atemp, b11, M2, half);

  // M3 = A11(B12 - B22)
  SubMatrix(b12, b22, btemp, half);
  Strassen(a11, btemp, M3, half);

  //M4 = A22(B21 - B11)
  SubMatrix(b21, b11, btemp, half);
  Strassen(a22, btemp, M4, half);

  //M5 = (A11 + A12)B22
  AddMatrix(a11,a12,atemp,half);
  Strassen(atemp, b22, M5, half);

  //M6 = (A21-A11)(B11+B12)
  SubMatrix(a21,a11,atemp,half);
  AddMatrix(b11,b12,btemp,half);
  Strassen(atemp,btemp,M6,half);

  //M7 = (A12-A22)(B21+B22)
  SubMatrix(a12,a22,atemp,half);
  AddMatrix(b21,b22,btemp,half);
  Strassen(atemp,btemp,M7,half);

  /* Now we have the Ms get the Cs */
  //c11 = M1 + M4 - M5 + M7;
  AddMatrix(M1,M4,atemp,half);
  SubMatrix(M7,M5,btemp,half);
  AddMatrix(atemp,btemp,c11,half);

  //c12 = M3 + M5
  AddMatrix(M3,M5,c12,half);

  //c21 = M2 + M4
  AddMatrix(M2,M4,c21,half);

  //c22 = M1 - M2 + M3 + M6
  SubMatrix(M1,M2,atemp,half);
  AddMatrix(M3,M6,btemp,half);
  AddMatrix(atemp,btemp,c22,half);

  /* Now fill in C */
  for(i = 0 ; i < half ; i++) {
    for(j = 0 ; j < half ; j++) {
      c[i*size + j] = c11[i*half + j];
      c[(i+half)*size + j] = c12[i*half + j];
      c[i*size + (j+half)] = c21[i*half + j];
      c[(i + half)*size + (j + half)] = c22[i*half + j];
    }
  }

  /* Free Everything */
  free(a11);
  free(a12);
  free(a21);
  free(a22); 
  free(b11);
  free(b12);
  free(b21);
  free(b22);
  free(c11);
  free(c12);
  free(c21);
  free(c22); 
  free(M1);
  free(M2);
  free(M3);
  free(M4);
  free(M5);
  free(M6);
  free(M7);
  free(atemp);
  free(btemp);

  return;
}

void printMatrix(float* a, int size) {
  int i =0, j =0;
  for (; i < size ; i++) {
    j = 0;
    for(; j < size ; j++) {
      printf("%f ",a[i*size + j]);
    }
    printf("\n");
  }
}
